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Mixed quantum-classical methods provide powerful algorithms for the simulation of 
quantum processes in large and complex systems. The forward-backward trajectory 
solution of the mixed quantum-classical Liouville equation in the mapping basis [J. 
Chem. Phys. 137, 22A507 (2012)] is one such scheme. It simulates the dynamics 
via the propagation of forward and backward trajectories of quantum coherent state 
variables, and the propagation of bath trajectories on a mean-field potential deter- 
mined jointly by the forward and backward trajectories. An analysis of the properties 
of this solution, numerical tests of its validity and an investigation of its utility for 
the study of nonadiabtic quantum processes are given. In addition, we present an 
extension of this approximate solution that allows one to systematically improve the 
results. This extension, termed the jump forward-backward trajectory solution, is 
analyzed and tested in detail and its various implementations are discussed. 



I. INTRODUCTION 



Nonadiabatic processes are important for the description of many chemical and biologi- 
cal processes such as proton and electron transfer reactions, vibrational relaxation, quantum 
reaction dynamics, photochemical dynamics, and coherent energy transfer phenomena in bi- 
ological systemd^. In all these processes, the dynamics of the system is strongly influenced by 
the environment in which it resides; hence, an accurate prediction of the system's properties 
entails the simulation of the system along with its environment. To avoid the exponential 
growth in the computational costs of exact quantum simulations with the number of de- 
grees of freedom, several mixed quantum-classicaP^ and semi-classical^ schemes have been 
developed. 

Simulation algorithms based on the quantum-classical Liouville equatiorP have been used 
to study complex open quantum systems that display nonadiabatic phenomena!^. Differ- 
ences among the algorithms often can be attributed to the basis set used to represent the 
quantum system, since the way the quantum degrees of freedom are propagated depends 
on basis set representations. Entangled trajectories^ can be used to simulate the quantum 
dynamics in the subsystem basis; surface-hopping-like algorithmd^^HI] g^n be formulated in 
the adiabatic basis; the force basi^^^*^, obtained from diagonalizing the Hellmann-Feynman 

force, yields yet another propagation scheme. Another class of simple and efficient al- 
gorithmpnSl 

can be constructed if the subsystem basis of an A^-level quantum system is 
mapped onto the single-excitation space of fictitious harmonic oscillators. In this map- 
ping representatiorP^^ the discrete quantum degrees of freedom are mapped onto continu- 
ous variables, such as the positions and momenta of the fictitious oscillators. The Poisson 
bracket mapping equationf^^'^^EQ] g^^^ ^-j^g forward-backward trajectory solutiorP^ are two 
such mapping-based approximate solution schemes for the quantum-classical Liouville equa- 
tion. The forward-backward trajectory solution, the topic of this paper, is derived from 
the forward-backward formpi' of the formal solution of the mapping-transformed quantum- 
classical Liouville equation. Since the quantum time evolution operator is expressed in a 
coherent state basis, the time evolution of a mixed quantum-classical system is simulated 
through the forward and backward trajectory evolution of coherent state coordinates 
while the bath coordinates evolve under the infiuence of a mean potential that depends on 
these forward and backward trajectories. In this scheme, the quantum-classical Liouville 
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dynamics is simulated via an ensemble of independent Newtonian trajectories. 

In this paper we analyze various aspects of the forward-backward trajectory solution. As 
shown earliei'^, this solution satisfies the differential form of the quantum-classical Liouville 
equation and is formally invariant to the form (trace versus traceless) of the Hamiltonian. 
To derive such a continuous trajectory picture, an orthogonality approximation is made, 
which assumes that the coherent states connecting various time segments are orthogonal 
to each other. Here, we show that the formal invariance to the form of the Hamiltonian 
is broken due to the use of the orthogonality approximation. If the system-bath coupling 
is not a weak perturbation to the pure bath potential, the use of the trace-form Hamilto- 
nian may yield poor results due to the possible presence of an inverted potential for the 
bath coordinates. To systematically improve the forward-backward trajectory solution, we 
generalize the simulation scheme by restricting the use of the orthogonality approximation 
to only certain time steps. In the resulting new solution, the forward and backward co- 
herent state trajectories experience discontinuous jumps in the phase space whenever the 
orthogonality approximation is not used; hence, the generalized solution is termed the jump 
forward-backward trajectory solution. Sampling of jumps is carried out by a Monte Carlo 
procedure. In order to improve the convergence, the focusing approximatiorp2H2l jg used in 
order to avoid performing the integrals over intermediate coherent state variables. 

The outline of the paper is as follows: In Sec. |Tl]we summarize the principal elements of the 
forward-backward trajectory solution and present its generalization that involves jumps. The 



focusing approximation is discussed in Sec. HI, and a simple picture of trajectory dynamics 



is given to illustrate the effects of focusing. In Sec. |IV| we present simulation results on 
four test models, which are chosen to highlight important aspects of the algorithms. The 
conclusions are given in Sec. |Vj 

II. FORWARD-BACKWARD TRAJECTORY SOLUTION 

We consider a quantum subsystem coupled to a bath where the dynamics is described by 
the quantum-classical Liouville equation (QCLE). The Hamiltonian has the form, 

Hw{X) = HbiX) + hs + VciR) = HhiX) + h{R), (1) 

where the subscript W refers to a partial Wigner transform over the bath degrees of free- 
dom (DOF). Here Hh{X) = P^/2M + Vb{R) is the bath Hamiltonian with Vb{R) the bath 



potential energy, hg = /2m + Vg is the subsystem Hamiltonian with p and Vg the subsys- 
tem momentum and potential energy operators, and Vc{R) is the coupling potential energy 
operator. The masses of the subsystem and bath particles are m and M, respectively. 

For a quantum operator which depends on the classical phase space variables 

X = {R, P) = R2, ...,Rni,, Pi, P2, -P/vJ of the bath, the formal solution of the QCLE 
is given b}^, 

Bw{X,t)=e'^'Bw{X), (2) 

where the QCL operator is 

i . . 1 . . . . 

iCBw = -[Hw,Bw] - -{{Hw, Bw} - {Bw,Hw}), (3) 



[Aw, Bw] = AwBw - BwAw is the commutator, and {Aw, Bw} = - ^ 
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the Poisson bracket with respect to X. We may also write the QCL operator a: 

zCBw = ^(^-Ha Bw - Bw Ha) , (4) 
which involves the forward and backward evolution operators, 

nA = Hwil + —j, nA=il + ^jHw, (5) 

with A the negative of the Poisson bracket operator, A =Vp ■ Vr — Vj? ■ Vp- In this case 
the formal solution of the QCLE can be expressed as 

Bw{X, t)=S (^e''^^'/^Bw{X)e~'^^'/^^ . (6) 

The S operator in this form simply specifies the order in which products of the left and 

right operators act in order to be identical with the first form in Eq. ([s]). In particular, a 

general term S [{l-iAy Bw{1-LaY^ in the expansion of the exponential operators is composed 

of (j + k)\/j\k\ separate terms each with a pre-factor oi j\k\/{j + A;)!. Each of these separate 

— > <— 

terms corresponds to a specific order in which the and 'Ha operators act on Bw- 



A. Hamiltonian dynamics in the coherent state phase space of the mapping 
representation 

We suppose that the time evolution of the quantum subsystem (coupled to the bath) can 
be accurately described within a truncated Hilbert space of dimension N . Furthermore, the 
subsystem basis {|A)} is chosen for the matrix representations of quantum operators. 



In the mapping representation, the state |A) is replaced by \mx), an eigenfunction of the 
Hamiltonian for fictitious harmonic oscillatord^'^, having occupation numbers which are 
limited to or 1: |A) Im^) = |0i,--- ,1a, ■■■Oat). Creation and annihilation operators 
on these states, dt^ and oa, satisfy the commutation relation [aA,dy] = 5a,a'- The actions of 
these operators on the single-excitation mapping states are d^ |0) = \rax) and a\ \m\) = |0), 
where |0) = |0i . . . Oat) is the ground state of the mapping basis. 

We may then define the mapping version of operators, Bm(X), given by Bm(X) = 
B^' {X)d\dy, such that matrix elements of Bw in the subsystem basis are equal to the 
matrix elements of the corresponding mapping operator: B^' (X) = {X\Bw(X)\\') = 
{■mx\Bm{X)\'m\i). (The Einstein summation convention will be used throughout although 
sometimes sums will be explicitly written if there is the possibility of confusion.) In partic- 
ular, the mapping Hamiltonian operator is 

H,n = H,{X) + h^^\R)d\dx, = H,{X) + h^, (7) 

where we applied the mapping transformation only on the part of the Hamiltonian that 
involves the subsystem DOF in Eq. The pure bath term, Hb{X) in Eq. acts as an 
identity operator in the subsystem basis and is mapped onto the identity operator of the 
mapping space. The formal solution in Eq. ^ may be expressed in mapping operators and 
is given by 

BUX, t)=S (e^«i'*/^5„(X)e-^«"A , (8) 

where "H™ is given by "H^ = Hm{^ + hA/2i), with an analogous definition for "H^. 

We now define the coherent states \z) in the mapping space, dx\z) = zx \ z) and {z\a\ = 
z^^ {z\, where \z) = \zi, . . . , zn) and the eigenvalue is zx = {qx + ipx)/V2^- The variables 
q = (gi, . . . , gTv) and p = (pi, . . . ,Pn) are mean coordinates and momenta of the harmonic 
oscillators in the state l^;), respectively; i.e., we have {z\ qx\z) = qx and {z\ px \z) = px- The 
coherent states form an overcomplete basis with the inner produce between any two such 
states, {z\ z') = _ Finally, we remark that the coherent states provide 

the resolution of identity, 

1 = / ^ 1^) (^1 ' (9) 
where d'^z = d{^{z))d{'^{z)) = dqdp/{2h)^. 

Following our earlier work?^, we then decompose the forward and backward evolution 
operators in Eq. (|8| into a concatenation of M short-time evolutions with Ati = t and 



Mr = t. In each short-time interval Atj, we introduce two sets of coherent states, \zi) and 
Iz'j), via Eq. (|9|. After some algebra, the matrix elements of Eq. ([s]) can be approximated 



by 



(mA {z[\ nix') 



MM 

e 



StV^'(X)(m^,|4f(AiA/))) (10) 
(^|...|4(At2)))(4kl(A^i)))- 



Equation (10) is evaluated sequentially, from smallest to largest times, by taking the bath 
phase space propagators in expressions such as e^^<^(^'^^'4)^'^i(^. ■ ■ ) to act on all quantities in 
the parentheses, including other propagators at later times, Atj>j. The bath phase space 
propagator reads iCeiX, z, z') = ^ ■ - '^^'^^"''^ -* ■ ^ with 

VeiX, z, z') = V,{R) -J2h'' + — {zlzy + ztz'y ) , 

A 

^ Vo{R) + I (h'^'zlzy + h'^'z'^z'y) , (11) 

= Vo{R) + ^ {KliR, z) + hd{R, z')) . 

In the first two lines, we suppress the /^-dependence of the matrix elements h'^'^' (R). Similar 
to other QCLE solutions, the bath dynamics are governed by classical trajectories in the 
bath phase space. 

To obtain Eq. ([To|, we also use the fact the coherent state evolution under a quadratic 
Hamiltonian hm{R) can be exactly evaluated by e"*'^™^ \z) = \z{Ati)), where the trajectory 



evolution of zx is governed by = —jr — . Thus the quantum subsystem dynamics 



dzx _ i dhci(R,z) 

is governed by segments of coherent state trajectories (defined by Zi and 4); which are not 



necessarily continuous across the time steps because the coherent state variables Zi and Zi+i 
are independent to each other. By making the orthogonality approximation to the inner 
products, {zi(ti) l-Zj+i) ~ n^S^Zi^i — Zi(ti)), we construct a smooth trajectory solution of 
the quantum dynamics. Furthermore, by scaling the coherent state variables 2;scaied = z/ V^, 



we also show that the trajectories in the scaled phase space (X, -Zscaied, -^^scaied) strictly 
Newtonian trajectories. In the rest of the article, we will consistently use scaled coherent 
state variables and drop the subscript "scaled". 

In summary, the orthogonality approximation and the scaling of coherent state variables 
helps to simplify Eq. (10) and its interpretation. We have 

B^'{X,t) = J2j dxdx'(l){x)(j){x') 

x^iq^{t)-tp^it))iq'^,{t)+tp'^,{t)), (12) 

where x = {q,p), dx = dqdp, and = (Z;,)"^ e"^"*^'^"^^^^/^ is the normalized Gaussian 
distribution function. We have written this equation in the scaled coherent state variables 
^\ = ilx + ipx)/^- The trajectories of Xt, Xt and x[ are governed by Hamilton's equations, 

dq^ _ dHe{X,x,x') dp^ _ dHe{X,x,x') 
dt dp^ ' dt dq^ 

dq'^ _dH,{X,x,x') dp'^ _ dH,{X,x,x') 



dt dp'^ dt dq'^ 

dR _ P dP _ dHe{X,x,x') 
~dt ~ M' ~dt ~ M ' 



(13) 



where (X, x,x') = P^ /2M + (R) + (^a^a' + PxPx' + qWx' + P'xP'x' ) • In the following 

discussion, we take h = 1 for simplicity. 

Similar to the trajectory solution of the Poisson bracket mapping equation (PBME), the 
FBTS approximates the quantum-classical dynamics in terms of an ensemble of independent 
Newtonian trajectories. However, the fact that the quantum-related phase space in the 
FBTS is twice as large as the PBME phase space allows more complex evolutions of these 



trajectories and more accurate characterization of the bath potential surface. In Sec. HI A 



we will discuss the reduction of FB trajectories to PBME-like trajectories under some specific 
conditions. 



B. Jump forward-backward trajectory solution 



As discussed in Ref. [15], Eq. (10) yields exact QCLE dynamics if one does not impose 



the orthogonality approximation on the coherent states. However, the large number of 
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intermediate variables Zi {z[) precludes practical computations of the integrals involving these 
variables. To extract QCLE results from FB trajectories, a systematic method is needed 
which can make an efficient compromise between computational costs and convergence. One 
possible method to achieve this compromise in a controlled manner is to select a subsequence 



of K time steps = 1,...,K} in Eq. (10) and evaluate the integrals of Zi^ and z'^^ 



explicitly {i.e. without resorting to the orthogonality approximations). According to this 
prescription, the continuous FB trajectories experience K discontinuous jumps in the (x, x') 
phase space. 

In practice, this scheme can be carried out in several ways depending on the selection of 
the K time steps. In the simplest case, one may select every (M/K) time steps from a total 
of M steps to fully evaluate the coherent state integrals in Eq. ((To|: 



K 

B^'{X,t) = J2 Yl Yldx,dx'J{x,)(f){x', 

UU' Sns'r,... v=0 



X (goA + tPoxMx' - tPox')B^''{Xt) (14) 



K 



X 



.v=l 
K 



X 



v=l 

X(gi^^(rx+l) - iPKf.i^K+l)){q'K^^'i^K+l) + ipK^,'{TK+l)), 



where the subscripts, v and s, refer to the f-th time step and the s-th component of q 
and p vectors, respectively, and = tj^ — with ti^ = and ti^^^ = t. To obtain 



Eq. (14), we inserted a projection operator V = Yls between the inner product of 

two coherent states at every selected time step, {zi{Ti+i)\V\zi+i) . In Appendix A, we show 
that the insertions of V do not introduce further approximations in this calculation. It is 
desirable to insert V every time we introduce a discontinuous jump to the trajectories, since 
then 0(xj) and appear at each time step. Otherwise, displaced Gaussian functions 

must be considered. The current form of the equation is also more suitable for introducing 
the focusing, an approximation scheme which is discussed in the next section. 

A stochastic generalization of the fixed-interval selection method can also be used. We 
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may write 



/ „ L 

UU' sns',,... v=0 



(15) 



X 



X 



V 

where = or 1. A prime (') is put on the second summation symbol and the last two 
product symbols in the equation to imply that the summations and multiplications only 
exist for a pair of and s'^ if k„ = 1. Furthermore, we replace the normalized Gaussian 
functions (p^x^) with 

S{xy — Xy^iir^)) , = 0, 



— 1. 



(16) 



The interpretation of Eqs. (15) and (16) is straightforward. We select L time steps that are 



J steps apart, i.e. JL = M. For a given binary sequence of . . . , kl}, we fully evaluate 
the coherent state integrals at the f J-th time step if = 1. For other selected time steps 
corresponding to k^, = 0, we still apply the orthogonality approximation to simplify the 
calculation. Finally, we remark that matrix element of interest, B^' {X,t), can be obtained 
as an average of all possible binary sequences of k, 



Bw i^^t) 



E 

ki,...kl 



Pi.}BZ..,.,{X,t\ 



(17) 



where P{k} denotes the discrete probability distribution of a given binary sequence of 



{ki, . . . ,kl}. Equation (17) reduces to Eq. (14) if one takes J = M/K and Pj^^^i} = 1. 



Also, Eq. (17) reduces to Eq. (12) if one takes P{k„=o} = 1 • 



We define the jump forward-backward trajectory solution (JFBTS) by Eqs. (15) - (17). 
In particular, we use the binomial distribution for Pj^} in this study. We take the Bernoulli 



trial probability p 



K 

M/J 



where K is the average number of jumps (equivalent to the average 



number of I's in the binary sequence of k). Equation (17) also has a simple interpretation in 
terms of trajectories if the integrals are evaluated using Monte Carlo (MC) sampling. During 



the propagation of a trajectory, a uniformly-distributed random number ^ is generated at 
every J-th step. If ^ < p, then we introduce a jump; otherwise, we continue to propagate 
the continuous trajectory until the next J-th step. 

Finally, we note the close resemblance in structure between the JFBTS and IPLDM ^ 
(iterative partial linearized density matrix), a similar generalization scheme that is applied 
to the PLDM algorithm in order to improve the results. The differences between JFBTS 
and IPLDM can mostly be attributed to the differences between FBTS and PLDA/PSIEZI 



III. FURTHER APPROXIMATION: FOCUSING 

In this section, we provide a detailed discussion of the effects of focusing, an approxi- 
mation proposed^^ by Bonella et ai, on the FBTS and JFBTS. Focusing can be defined as 
follows, 

dx(f){x){ql+pl)f{x) 

dx6iql+pl-l)ll5iq,)6ip,) 

xiql + pDfix), (18) 

where f{x) should be at most weakly dependent on x variables. In order to obtain the 
approximate integral, one defines the mapping weight p(x) = (g^ +p\)(f){x), and assumes 
that only the maximum points, satisfying the Dirac delta functions in Eq. ( jlsj ), of this weight 
contribute significantly to the integral. 



Equation (18) may not seem to be directly related to any integrals we have presented 



so far. However, if one considers the integral in Eq. (12) with t ^ 0, then the integrand 
'^/ifj.' i.-^t){,qii{t) — Wii{i)){.Q'pL'{'^) + can be approximated by a single term (where 

fi = X and ^' = A') of the summation. Furthermore, near the initial time t = 0, the 
integrand B^' (Xt) is weakly dependent on coherent state variables, and the polynomial 
factors involving x and x' can be approximated by their initial values. Same arguments are 



also applicable to Eq. (15). Hence, focusing can approximate the integrals of Eqs. (12) and 
(15) fairly well when t is small or the system-bath coupling is weak; i.e., a separation of 
time scale can be assumed. Unfortunately, these assumptions are not valid in the long-time 
limit for most physical systems. 
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By using focusing, one avoids the full sampling from the Gaussians in the MC evaluation 
of Eqs. (12) and (15). As indicated in many numerical test j-^^ * ^^ * ^^ * ^^ !, focusing often yields a 
converged result with at least an order of magnitude fewer trajectories than that required for 
a full calculations. However, when used indiscriminately, focusing might yield poor result^. 



A. Focused Initial Condition in the FBTS 

We suppose that the quantum subsystem is initially in a single subsystem state |A). This 
requirement is not as restrictive as it might appear, since our formalism does not dictate 
which particular quantum subsystem basis is to be used. Therefore, we can often take the 
initial state of the system to be a vector in the basis. When focusin^^S jg applied to the 
FBTS, the approximation restricts the initial conditions of the trajectories to a very narrow 



region in the phase space specified by the Dirac delta functions in Eq. (18). 

First, we discuss the effects of the focused initial conditions on the forward and back- 
ward coherent state trajectories. In Appendix B, we prove that the backward coherent state 
trajectories can be replaced by the forward trajectories in the FBTS if focused initial condi- 
tions are imposed and the required assumptions (such as the separation of time scales) are 
valid. This observation shows that the two sets of trajectories are highly synchronized in the 
sense explained in Appendix B. One manifestation of such synchronization is the reduction 
of noise in the calculation of the phase factor of the quantum subsystem density matrix 
element, p^^'(X, t). In the full FBST calculation, such a matrix element is determined by 
the ensemble average of {qx + ip\){q'x — — Wti{t)){.l'^i{i) + w'^'if)) according to 

Eq. ( [I2| ). However, in the case of focused initial conditions, one can simply replace {q',p') 
with {q,p) in the above expression for the ensemble average. For instance, if one employs 
the formalism to compute a diagonal density matrix element, p^^(X, t), then one finds that 
this replacement yields an exact real number as required. For the full calculation, one needs 
to average over a large number of trajectories to suppress the noisy imaginary components. 

Next we consider the effects of focused initial conditions on bath trajectories. As men- 
tioned earlier, focusing should be a valid approximation for an initial short time interval or 
for a weak subsystem-bath coupling. We will take this into account by assuming that spatial 
derivatives of C{R), which is a column matrix of eigenvectors of h{R), can be neglected. 
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Under this assumption, the bath momenta evolve according to, 
dP dVo{R) dh^^\R) 



dt dR dR 



iqu{t)qu'{t)+p,{t)p,it)) 



X cos(Aa;^^/t) (g^ + p2 j 

X cos(Aw^^/t), (19) 

where Au^^' = — uj^i. In writing this equation, as discussed above, we replaced the 
backward coherent state variables with the forward ones. To obtain the second equality, 
the equation of motion for the coherent state variables was integrated to yield analytical 
expressions for q{t) and p{t). The final form of this equation was obtained using the focused 
initial condition p\ = In the full calculation, the system-bath coupling term (the 
second equality in the equation) should be averaged over all possible values of X]x(^x~'~^x) ~ 
r with respect to the Gaussian distribution (j){r). 

For situations typical of a scattering problem, where the quantum subsystem is initialized 
in the asymptotic region with C^y = (i.e. the subsystem basis and the adiabatic basis 



coincide), the final form of Eq. (|l9j) reduces to dP/dt = -dVoiR)/dR - dh^^{R)/dR. This 
simple bath dynamics is the adiabatic approximation in the asymptotic region. In contrast, 
the full calculation requires summation over a large number of trajectories such that the 
interferences among them leads to the adiabatic approximation. However, in the interaction 
region of the scattering potential, this adiabatic approximation is certainly not sufficient. 



Even though Eq. (19) shows that focusing includes more than the adiabatic approximation 
on the bath dynamics, it simply does not include enough corrections to fully describe the 
dynamics in the interaction region in many cases. For instance, in a recent calculations on a 
simple avoided crossing modeP^, Huo and Coker showecP^ that the correct bath properties 
cannot be reliably extracted from focused initial conditions. 



B. Iterative Focusing in the JFBTS 



While focusing is not an essential part of the FBTS, it is often the only tractable means by 
which one can obtain a converged JFBTS solution when one needs more than five jumps in 
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order to obtain the QCLE result. Similar to several iterative schemes of partially linearized 
path integral formalism^^SEl^ such as IPLDM, we rely on focusing's superior convergence 
properties to minimize the increasing number of trajectories required to produce a stable 
result. 



Comparing Eqs. (12) and (15), the JFBTS essentially adds a sequence of summations over 



intermediate quantum states, Sj and s^, and integrals over coherent state variables Zi and z[ 



on top of Eq. (12). While an exact evaluation of Eq. (15) is conceptually straightforward, 
the numerical implementation incorporating the focusing approximation is slightly more 
complicated. First we evaluate the sequence of intermediate summations by an importance 
sampling MC calculation. This technique has been adopted in several related QCLE and 
linearized path integral algorithm j^^ * ^^ * ^^ l Thus, in each MC run, we randomly select a pair 
of forward and backward quantum states (sj, at every i-th jump instead of summing over 
all quantum states. The pair selection is governed by a discrete probability distribution, 

Ps.,s'^ = {ql + vDiq'l + Kp/M with M = Y^s^^qI + vDiq'l + Kp- This probability 
distribution is based on the intuition that the selection of the pair indices should be related 
to the present distribution of the population weight, (g^ over all quantum states. A 

re-weighting factor, l/Pg- is added to the MC evaluation whenever such an importance 
sampling is made. Once the selection of a pair of quantum indices is made for the i-th 
jump, we use the focusing approximation with the condition + — ^)5{q'^i + v's' ~ 1) 
to evaluate the integral over Zi and z[. In the above focusing condition, we imply that all 
unspecified components of forward and backward variables are set to zero. 

Now we analyze the effects of iterative focusing on the trajectories in the JFBTS formal- 
ism. First, we provide an intuitive picture of the trajectory evolution when full sampling 
over coherent state variables and full summation over quantum states at each intermediate 
step are performed exactly. At the initial time, t = 0, we assign initial conditions of 
and Pi, for each quantum state |mj,) by sampling from the Gaussian The trajectories 



are then propagated with Hamiltonian dynamics based on Eq. (13). At each time step, we 
stochastically decide if a jump should be introduced. When a jump is introduced at time t, 
we store a snapshot of the trajectories by saving the polynomial factor {qv{t) ± ipv{t)) for 
each z/, where the plus and minus signs are for forward and backward coherent state vari- 
ables, respectively. Then we discard the trajectories and re-sample new initial conditions for 
those trajectories. We also take a snapshot of the new trajectories by saving the polynomial 
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factor (qi, =f ip,^) prior to further propagation. The same process is then be repeated until 
the simulated time length has been reached. Due to the interferences in the full summations 



of the products of the polynomial factors in Eq. (15), only those trajectories with rather 
minimal discontinuities contribute significantly to the final result. 

Next we describe a rather different evolution of trajectories in the phase space. The 
combination of the importance sampling of intermediate quantum state indices and the 
focusing approximation for the coherent state integrals yields a picture involving trajectory 
collapse whenever a jump is introduced. In this numerical implementation, each forward and 
backward trajectory can only be initialized in the given pair of quantum states, \ms-) and 
|ms^). Once the initial conditions are assigned to the trajectories, the population weights, 
Qs + Ps) of the trajectories will gradually shift from the initial states, Ims^) and to 
other states in the forward and backward phase space, respectively. After some time, the 
population weights may become more uniformly distributed over all available states. There 
will then come a moment where a jump is introduced, a new pair of quantum indices is 
selected, and a new focusing condition is imposed to collapse all trajectories onto the new 
pair of forward and backward states. 

These two rather different pictures of the trajectory evolution make one wonder if iter- 
ative focusing applied to the JFBTS can be trusted. We argue that the result should be 
reliable. Different from the focused initial condition where the states to be focused onto 
are fixed, iterative focusing collapses the forward and backward trajectories into all possible 
combinations of pairs of quantum states. Therefore, if we have a large enough MC sample, 
we will eventually obtain collapsed trajectories re-initialized in all possible quantum states 
and recover an ensemble picture that is in accord with that of the full-sampling case. 



IV. SIMULATION RESULTS 

We now present the results obtained from simulations of the FBTS and JFBTS for four 
models, which are selected to examine different aspects of the simulation algorithms. In 
particular, we use these results to illustrate the following: First, the differential form of the 
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FBTS has been shown!^ to be invariant to the following forms of Hamiltonian, 



Hw(x) = < 



Hb{X) + h{R), trace form, 

(20) 

Hb{X) + h{R), traceless form, 



where Hb{X) = Hb{X) + Tr,/i/iV and h{R) = h{R) - TTsh{R)/N where Tr, denotes the 
trace over subsystem DOF. The local, infinitesimal time-step analysis can be extended to 
a finite-interval time scale if the orthogonality approximation is not used and, in this case, 
results using either form of the Hamiltonian will be identical. However, if the orthogonality 
approximation is made, global errors which depend the form of the Hamiltonian will accu- 
mulate. The traceless version of the Hamiltonian often yields similar or better results than 
the trace form, since it is less susceptible to problems due to inverted potentials. Analysis 
of the results of the dual avoided crossing (Tully 2) and the Fenna-Mathews-Olsen (FMO) 
models support this conclusion. For consistency, we will display the Hamiltonian in trace- 
form when we introduce the models in this section, although calculations have been carried 
using both forms as indicated below. 

Second, we analyze the effects of focused initial conditions on the FBTS. We use results 
of the Tully 2 model and spin-boson models to illustrate the trajectory picture described 



in Sec. Ill Through these examples, we show that the Tully 2 and symmetric spin-boson 
models represent a small set of special cases where the use of focused initial conditions might 
yield results similar to those of the full sampling calculation. 

Last, we consider the application of iterative focusing with JFBTS to alleviate the onerous 
demands on the very large number of trajectories required in the full-sampling calculations. 
From the results on the asymmetric spin-boson and conical intersection models, we demon- 
strate the convergence of the JFBTS to the QCLE solution. Due to the fact that trajectories 
collapse onto a pair of diabatic surfaces at each jump, it turns out that the selection of the 
subsystem basis {i.e., which diabatic surfaces are used) in the formalism affects the conver- 
gence of the results. We provide general guidance on how to select an ideal subsystem basis 
for a JFBTS calculation. 
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A. Dual Avoided Crossing Model 



The subsystem Hamiltonian for the Tully 2 model in the diabatic basis {|1) , |2)} is 

Ce-^^' 

(21) 

The numerical values of parameters A, B, C, D, Eq and all other details of this particular 
model calculation are taken from Ref. [I9]. The partially Wigner transformed Hamiltonian 
of this system is Hw{X) = P^/2M + h{R). Note that the pure bath term contains only the 
kinetic energy of a single bath particle. Initially, the quantum subsystem is in the state |1) 
and the bath particle is modeled as a Gaussian wave packet centered at Rq with initial bath 
momentum Pq directed towards the interaction region. 

Figure [l] shows the asymptotic populations of the quantum subsystem in the two diabatic 
states as a function of the initial momentum Pq. In this figure, we compare three different 
FBTS results (traceless form, trace form, trace form with focused initial conditions) with 
the exact quantum results. At high momentum, Pq > 35, all different FBTS results converge 
to the exact quantum calculation. This convergence is due to the fact that Gaussian wave 
packet passes through the interaction region of the scattering potential with high velocity, 
and the system-bath coupling does not strongly influence the quantum and bath dynamics. 
In the low momentum regime, the system-bath coupling plays a much more crucial role in 
the bath dynamics, which in turn influences the quantum dynamics. It is clear that the 
traceless form (light blue curves) result match to the exact result the most. 

We expect that simulations using the traceless-form of the Hamiltonian will yield the 
best results, due to the possible presence of an inverted potential when the trace-form 
Hamiltonian is used. The fact that the Tully 2 model has no pure bath potential exacerbates 
the errors arising from dynamics on an inverted potential surface. It is interesting that 
the trace-form result improves significantly when the focusing approximation is used. We 
attribute the success of the focusing approximation in this case to the validity of the adiabatic 
approximation in the low momentum regime (up to roughly Pq = 15). The probability of 
finding the quantum subsystem in the adiabatic ground state at any time is over 90% 
for initial momenta up to Pq = 15. Since focused initial conditions describe adiabatic 
dynamics and some additional improvements, the much more accurate low-momentum result 
(in comparison to the full sampling calculation with the trace-form Hamiltonian) can be 
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h{R) = 



understood. 

Next we discuss how the FBST's formal invariance to the trace of the Hamiltonian is 
recovered as we introduce discontinuous jumps in the phase space. We studied the con- 
vergence of the asymptotic population as a function of the number of jumps introduced in 
the JFBST calculations with both forms of the Hamiltonian. As expected, we found that 
both sets (traceless and trace forms) of results improved and converged to each other. For 
instance, for the case of Pq = 20, the trace-form results converge to the exact result shown 
m Fig. [1] with 4 jumps. Similarly, the traceless-form results also converge to the exact result 
with 2 jumps. 

From these results, we see that the traceless-form FBTS is the preferred algorithm for 
numerical computations. We also see the usefulness of focused initial conditions in situations 
where the dynamics is weakly nonadiabatic. 

B. Spin-Boson Models 

The symmetric (unbiased) and asymmetric (biased) spin-boson models provide additional 
insight into the utility of the algorithms. Since spin-boson models and their generalizations 
have been studied extensively in different sub-fields of physics and chemistry, we focus on 
an explanation of the performances of FBTS and JFBTS algorithms in the reproduction of 
exact quantum results. 

The partially Wigner transformed Hamiltonian of the spin-boson model reads, 

+ea, - na^, (22) 

where Mj and Ui are the mass and frequency of the i-th bath oscillator, respectively, q 
controls the bilinear coupling strength between the i-th oscillator and the two-level quantum 
subsystem, Q is the coupling strength between the two quantum levels, e is the bias, and 
d'z(x) is the z{x) Pauli matrix. We assume that the bilinear coupling in the spin-boson 
models is characterized the ohmic spectral density, J(w) = vr c^/(2Mju;i)5(u; — Wj), where 
d = {^AMj^/'^Ui, Ui = -Wcln(l - iAu/uc), and Au = ujc{l - e-'^'""-/^-)/A^B with Uc the 
cut-off frequency, Nb the number of bath oscillators, and ^ the Kondo parameter. In all the 
spin-boson models considered here, the quantum subsystem is initialized in state |1) and the 
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bath is initially in thermal equilibrium. 

First we consider FBTS results for the symmetric spin-boson model. Figure [2] shows the 
time evolution of the population difference, {(Tz)- We not only consider the full-sampling and 
focused initial condition FBTS results, but also consider two cases with modified focused 
initial conditions. In case 1, we replace the typical focusing condition 5{q\ + — 1) by 
5{q\ +p\ — 1.5). In case 2, we use the even more severely modified focused initial condition 
6{ql +pI- 0.8)5(g| + pi - 0.2). In both modified cases, we omit the backward coherent 



states because they can be exactly replaced by the forward states as discussed in Sec. |III A 
All FBTS results in Fig. [2] are scaled in the following way: < cr^ > /lj\/c, where Imc is the 
MC estimate of the constant 1 using either the full sampling, proper or modified focused 
initial conditions, depending on the case considered. 

Examination of Fig. [2] reveals a rather important reason behind the reported account^^^'^ 
of high accuracy in the application of focused initial conditions to the symmetric spin-boson 
model. The results in the figure show that any modified focused initial condition of the form 
5{q\+p\ — r) with r being an arbitrary real number can yield results (scaled by Imc) almost 
identical to the exact quantum calculation. In addition, case 2, where the modified initial 
condition has non-zero quantum amplitude in state |2), also yields good results, except that 
the initial amplitude is slightly lower. This behavior, in which a large portion of phase space 
leads to similar quantum evolution, can be traced to the fact that the two diabatic states 
are energetically degenerate in the absence of the bath. 

For FBTS and related mapping algorithms, it is crucial to simulate the correct bath 
dynamics. We argue intuitively that the energetic degeneracy of the two states gives an 
effectively trivial quantum feedback on the bath. In the ensemble picture, the bath oscillators 
oscillate around their equilibrium points as if the quantum subsystem only influences the 
amplitudes of their oscillations but does not modify their oscillation frequencies or displace 
them from their equilibrium points. However, the quantum subsystem still feels the bath 
dynamics and suffers decoherence. We support this intuitive argument by the following 
observations: First, if one simply suppresses the subsystem-related part in the equation of 



motion for P in Eq. (13), i.e. use dP/dt = —dVo{R) /OR, then one still obtains a result which 
is almost identical to that shown in Fig. [2] Second, after the fifth oscillation in the figure, 
even case 2 converges to the exact result. This suggests that once the system moves beyond 
its correlation time scale with respect to the initial condition, the system is insensitive to 
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the population distribution over the two energetically degenerate levels. 

The asymmetric spin-boson model in which a bias e is introduced is a more challenging 
system for simulation. The consequences of having a bias in the spin-boson model have 
been analyzed extensively; for instance, the fundamental differences between the quantum 
dynamics in both symmetric and asymmetric spin-boson models have been summarized in 
Ref. [30]. 

Figure [3^ presents the time evolution of < cr^ > for several situations: FBTS (full sam- 
pling), FBTS (focused initial condition) and JFBTS (with 26 jumps). For the asymmetric 
spin-boson model, we find that the focused initial condition result deviates significantly 
from both the full-sampling and exact quantum results, after about Qt > 1.5. The rather 
poor performance of the focusing approximation to the FBTS can again be understood 
by analyzing the trajectories with modified focused initial conditions. In Fig. [3]d, sev- 
eral such cases are presented. For cases 1-3, the modified initial conditions take the form 
6{ql +pI — r) with r = 1.5, 1.8, and 2.0, respectively. In case 4, we use the initial condition 
6{ql + pf — 0.8)(5(g| + pi — 0.2). All the results in Fig. [3]d are scaled by the corresponding 
Imc as explained earlier. When e 7^ 0, different initial conditions yield very different results 
for ((Tz). Comparing Figs. [3^ and b, it is clear that the proper focused initial condition is 
simply not a valid approximation since it prunes too many trajectories from the ensemble. 

In Fig. |3^, we also present results for a JFBTS simulation with 26 jumps. We remark 
that the minimum number of average jumps required for a JFBTS calculation to recover the 
QCLE result depends on the interplay of several factors, such as the size of the time-step. 



the probability distribution Pj^} and the time-block size J discussed in Sec. II B In this 
calculation, we used the binomial distribution for P{k}, and a jump could be introduced in 
every even-number time step. Finally, we remark that, in certain cases, the QCLE result 
might be obtained with fewer jumps if one uses fixed-interval jumps. In summary, the 
JFBTS formalism is flexible in the sense that many parameters and the jump process can 
be modified in order to optimize its performance for each individual problem. 

From the analysis of these spin-boson models, we see that the applicability of focused 
initial conditions is very limited in scope. Their ability to describe symmetric spin-boson 
cases should be taken simply as a rare exception. However, the use of iterative focusing in 
the JFBTS was shown to be a more robust and useful tool. It allows us to obtain the QCLE 
solution, which agrees with the exact quantum result. 
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C. Fenna-Matthews-Olson Model 



The Fenna-Matthews-Olson (FMO) protein of green sulfur bacteria plays a role in the 
transfer of excitation energy to the reaction center involved in photosynthesis^. The model 
Hamiltonian for this process provides an example of a multi-level quantum system coupled 
to a bath. The Hamiltonian is composed of a seven-level quantum subsystem corresponding 
to an excitation state localized in one of the seven pigment proteins of the FMO complex, 
and each quantum level is bilinearly coupled to its own set of bath oscillators; the bilinear 
system-bath coupling is characterized by the Debye spectral density. The explicit expression 
for the Hamiltonian and the parameters used in its definition will not be given here but is 
available in the Ref. [32] and its supporting documents online. The quantum subsystem is 
initialized in the state |1) and all the bath oscillators are initially in the thermal equilibrium 
state. This model has been studied often. In particular, the PBMEp^ and the PLDM!^ 
algorithms, closely related to the FBTS, have already been applied to investigate the pop- 
ulation dynamics of the FMO complex and compared to the numerically accurate quantum 
result^™. 

Figure [i] plots the populations in states |2), and |3) as functions of time. Two sets of 
results are presented and compared against quantum results obtained from rescaled Hierar- 
chical Coupled Reduced Master Equation (HCRME) calculatiorP^I. The solid lines represent 
the results of the traceless-form FBTS, while the corresponding colored dots correspond to 
results of the trace-form FBTS. For this system the FBTS is rather insensitive to whether 
the trace or traceless forms of the Hamiltonian are used. Unlike the Tully 2 model, the FMO 
model has a dominant pure bath quadratic potential for each oscillator. Since the system- 
bath coupling is a weak perturbation, it is less likely to encounter an inverted potential for 
the bath trajectories. In addition to the agreement between the trace-form and traceless- 
form, the results are also in quantitative accord with the rescaled HCRME calculation^. 
Furthermore, we have also checked and found (not shown) that the long-time population 
distributions of these states approach the thermal equilibrium distribution at time t = 10 
ps. 

The results of this section demonstrate the efficiency of the FBTS in dealing with multi- 
level systems. For this system, inverted potentials in the simulation of trajectories do not 
present a problem and the FBTS is not sensitive to the form of the Hamiltonian. 

20 



D. Conical intersection model 



Finally we consider a two-level, two-mode quantum model for the coupled vibronic states 
of a linear ABA triatomic molecule constructed by Ferretti, Lami and Villiani (FLV). In 
this model, the nuclei are described by two vibrational degrees of freedom, X and Y. The 
partially Wigner transformed Hamiltonian reads 

(23) 

where the subsystem Hamiltonian is defined by the following matrix elements: 

h'\R,) = ^MxcoliX - X,f + A, 

h'^{Rs) = ^Yexp (-a(X - X^f - . (24) 

In these equations, Rs = {X, Y), Ps = {Px, Py), Mx,y and ojx,y are the mass and frequency 
for the X and Y DOF, respectively. The quantum subsystem is initialized in the adiabatic 
ground state, while the vibronic X and Y initial states are taken to be Gaussian wave 
packets. Further details of this model can be found in Refs. [33] and [55] . 

Figure |5] plots the adiabatic ground state population at t = 50 fs as a function of the 
coupling strength 7. We see that the FBTS result matches the Trotter-based QCLE (with 
filtering) result up to about 7 = 0.02 but deviates significantly from the QCLE and ex- 
act quantum results for larger coupling strengths. At higher coupling strengths, the errors 
induced by the orthogonality approximation become significant. The result of a 15-jump 
JFBTS calculation is also shown. This simulation is able to reproduce all trends in popula- 
tion versus coupling strength curve and, for this number of jumps, about 10 % accuracy is 
achieved at the highest coupling. 

Next we discuss the importance of the subsystem basis used for the JFBTS calculation. 
The JFBTS result in Fig. [sjwas obtained using a representation with the basis states |±) = 
(|1) ± |2)) /\/2 where subsystem states |1) and |2) were used to define the Hamiltonian 



matrix elements in Eq. (24). In the figure, this rotated basis yields a JFBTS result with 
as much as ~ 25% improvement over the FBTS results. In our other JFBTS calculation 
(not shown in the figure) with respect to the original basis, {|1) , |2)}, we obtained at most 
~ 15% improvement over the FBTS results with the same number of jumps. 
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As discussed in Sec. |III B[ the importance sampling used to select a pair of forward and 
backward state indices, for the collapse of the trajectories in each MC run helps 



to avoid doing the exact summation over state indices in Eq. (15). Thus, another key 
to the efficiency of the JFBTS algorithm is to have a dispersed probability distribution 
for Psi,s'. in order to select as many different combinations of states as possible. An easy 
method to satisfy this requirement is to simply use the uniform probability distribution 
Psi,s'- = where is the number of quantum states. However, we prefer the distribution 



Psi,s\ = (q's- +PsJ('?s^ +p's')/-M defined in Sec. IIIB since it minimizes the abrupt changes 
that can influence the bath dynamics. Another way to obtain a more uniform distribution is 
to select a basis where the off-diagonal matrix elements are maximized since these influence 
population transfer among different levels. 

In the original basis, the FLV model has a sharply peaked off-diagonal matrix element, 
h^'^{Rs). Thus, in the asymptotic region where h^'^^Rg) ^ 0, the population distribution 
over the quantum states does not change signiflcantly even if jumps are introduced in the 
JFBTS calculation. Figure [6^ shows the population distribution over pairs of forward and 
backward quantum states for a typical trajectory in the original basis. The distribution is 
highly localized on one particular pair of states, except when the system moves through 
the interaction region where h^'^{Rs) approaches maximum strength; however, h^'^{Rs) has 
a rather steep proflle and the interaction region occupies a very small portion of the phase 
space. Therefore, extensive MC sampling involving many jumps would be required to sample 
all different pairs of states. By contrast, in the rotated basis K^s) = ±/i^^ are 

now the diagonal elements of the Hamiltonian matrix, and the off-diagonal matrix element 
h~^~{Rs) is linearly dependent on X. In the "asymptotic" region (where h^'^{Rs) ^ 0), 
the two diagonal matrix elements are approximately energetically degenerate in the rotated 
basis. Consequently, population transfer occurs easily in the "asymptotic" region. Since the 
"asymptotic" region represents the bulk of the phase space. Fig. ^jp shows a more balanced 
and oscillating population distribution over pairs of forward and backward quantum states 
for a typical trajectory in the rotated basis. 

These results show that the FBTS formalism can be customized and several modiflcations, 
such as the choice of subsystem basis and the jump probability distribution, etc., can be 
independently tuned and adapted to suit the problem under study without actuaUy changing 
the structure of the algorithm. From the analysis of FLV model, we have seen how to 
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determine the most suitable subsystem basis for the JFBTS calculation. 

V. SUMMARY AND CONCLUSION 

In this work, we demonstrated the usefulness of the FBTS through simulations of a va- 
riety of systems. In most instances, we obtained converged results with ~ 10^ trajectories. 
These numerical tests suggest that the FBTS performs at the same level of accuracy and ef- 
ficiency as the PLDA/P^EMZ]^ closely related algorithm. We showed that the traceless-form 
Hamiltonian yields more accurate results when there is no dominant pure bath potential in 
the model and should be used for all numerical computations. In cases where the FBTS 
results were compromised by the orthogonally approximation, we demonstrated that these 
results could be systematically improved by increasing the number of jumps in the general- 
ized JFBTS algorithm. In order to make JFBTS a computationally efficient algorithm, we 
carried out a detailed investigation of the focusing approximation and provided an extensive 
analysis of the approximate trajectory dynamics. This analysis indicated that focused ini- 
tial conditions should be used only when the nonadiabatic effects can be treated as a weak 
perturbation or for special cases such as the symmetric spin-boson model where the bath is 
rather insensitive to the differences between the quantum states. Finally, we also justified 
the validity of iterative focusing in the JFBTS. In all cases studied the QCLE solution either 
reproduced the exact quantum results or was a very good approximation to the exact results. 
Thus, the ability of the FBTS and JFBTS to reproduce the QCLE solution also implies an 
ability to reproduce the full quantum results. 

The number of jumps used in the JFBTS determines how many trajectories will be 
needed to obtain a converged result. As noted throughout the paper, several implementation 
details of the JFBTS algorithm can be customized to significantly reduce the minimum 
number of jumps and improve the convergence properties of the algorithm. In addition, 
the number of jumps needed to obtain QCLE result will certainly depend on the details of 
the Hamiltonian of each model. It is not a simple manner to predict the number of jumps 
and trajectories required to obtain the full QCLE result for an arbitrary system. However, 
at least two orders of magnitude increase in the number of trajectories (in comparison to 
the corresponding FBTS calculation) should be expected for most cases in which the FBTS 
results are compromised by the orthogonality approximation. Perhaps the greatest utility 
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of the JFBTS algorithm is that it provides a simple method to gauge the sufficiency of the 
FBTS when one studies nonadiabatic dynamics for complex real systems, beyond the simple 
benchmark models studied here for which exact quantum results are known. 

Potential improvements to the family of (J) FBTS algorithms are still possible; especially, 
it might be helpful to reformulate the entire formalism in the adiabatic basis. Currently, 
the flexibility of the algorithm can be traced back to the use of the subsystem basis, which 
is not uniquely defined as in the case of adiabatic basis. In particular, a current challenge 
of the JFBTS algorithm is related to the determination of the best timings for the insertion 
of jumps in the course of trajectory evolution. In the current formulation, there are no 
indicators which help to make this judgement. However, as is often the case with surface 
hopping algorithms formulated in the adiabatic basis, we expect that inherent indicators 
naturally emerge from the formalism once it is re-expressed in the adiabatic basis. 
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Appendix A: Insertions of Projection Operators 

In this appendix we support the claim that insertions of projection operators, V — 
^ \ms) in the mapping representation do not introduce additional errors. We first 
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consider the concatenation of parameterized time evolution operators in the subsystem: 

X ...(/iM|e-sM^M-i)-|A') 



X . . . {nif.Je-^'^^'^^^'-^^^lmx'). (Al) 

On the right side of the first equahty, we insert the resolution of identity X = at 
every time step. To obtain the second equality, we use the matrix identity: ei'^^^^'^ = 
(m^l es''™("^)'^ Im^'). If we express every time evolution operator in terms of Zi variables, 

et^'"(^)^= / ^k)(^(r)|, (A2) 



where {z{t)\ is time-evolved by Eq. (13), we obtain the desired expressions in Eq. (14). 



Appendix B: Forward and Backward Trajectories with Focused Initial 
Conditions 

In this appendix we take h = 1 and use the scaled coherent state variables as defined 
in Sec. [Tlj We prove that the backward trajectories can be exactly replaced by the forward 
trajectories in the FBTS if the focused initial condition is imposed. The coherent state 
variables, z and z', matter critically in two places in the formalism. First, the equation 



of motion for the bath momenta defined in Eq. (13) depends explicitly on these variables 



through the term | gjj^^^ {qxQx' + PxPv + qWx' ~^ p'xP'x')- Second, these variables play a 
critical role in the evaluation of the expectation value of an operator B, 



{Bit)) = Y.h^ I dxdx'B^\X,t)p^'iX), 

X,X' 

= Y [ dX f dxdx'^{x)(p{x')p^'{X) 



X,X' 

x{qx-^Px){q^.{t) + ^p^{t)) (Bl) 
xiq'y + tp',,){q'^,{t)-tp'^,{t))B!^^'{X,l 
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where the q and p variables with no exphcit time arguments are taken to be the initial values 
at t = 0. In the following, we prove the identities: 



<l'x{tyx'it)+Vxit)Mt) = qx{t)qx'{t)+px{t)py{t), 

{q'x' + w'x'W^'i^) - W) = (^v + ipx'){qt,'{t) - iPf,'{t)). 



(B2) 



The first identity implies that the backward trajectories can be replaced by the forward ones 
in the equation of motion for the bath. The second identity implies that the expectation 
value of the operator B{t) can also be evaluated without having the backward trajectories. 

Without loss of generality, we assume the quantum state is initialized in state |1) with the 
following focused initial condition: 6{ql +pl — l)6{q'i +p'i — 1) and every other component 
of the q and p vectors is set to zero at t = 0. We remind that the q{t) and p{t) variables 



satisfy the equation of motion in Eq. (13). First we define a function gxx'{t) = (l'x{t)Q'x'i^) + 
P'xi't)p'x'{t) - qx(t)qx'it) - Px{t)px'{t). It is easy to verify that gxx'{^) = for every (A, A') 



combination. Now we prove that 



d"9xx'it) 



for all n. We show that the first two 



derivatives of gxx'(t) are equal to zero: 
dgxx'{t) 



t=o 



dt 



h ""[pWx' - QaP'x' -PaQx' + qaPx']\t=o 

+h''''[pWo 
0, 



QxPa 



Pxqa + qxPa\\t=0 



(B3) 



and 



d'gxx'it) 



2h''^h''^g^f,{0) - h'^'^h^^gpxiQ) - h>^^K^^ gp^,{Q) 



t=o 



0. 



(B4) 



Since any higher derivative of gxx'if) at t = can be recursively defined in terms of Eqs. (B3) 



and (B4), we have proved that gxx'{t) = for all t and established the first identity in 



Eq. (B2) 



Next we define functions fxx'{t) = q'xq'x'{i) + P'xP'x'i^) - QxQx'ii) - PxPx'{t) and kxx'{t) 



~q'xP'x'{^) ~^ PWx'it) ~ Pxqx'{t) + qxPx'{t). We remark that the second identity in Eq. (B2) is 
equivalent to fxx'if) +ik\\i{t) = 0. Similar to the previous analysis, we note that fxx'it) = 
and kxx>{t) = by showing that these functions and all orders of their derivatives are zero 
at t = 0. Furthermore, we note the simple relations between the first derivatives of these 



26 



functions: 



dfxx'it) 



dt 



dkxx'(t) 



h^'''WxP'a{t)-pWa{t) - qxPa{t)+Pxqamt=0, 



t=0 



0, 



dt 



t=0 



h^''[qWa{t)+PxPa{t)-PxPa{t) " gAg«(t)]|o, 

/^''"/a.(0) 
0. 



(B5) 



From these equations, it is obvious that all higher derivatives of these functions evaluated 



at t = should vanish. Thus, we have established the second identity in Eq. (B2) 
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FIG. 1. Asymptotic populations of the diabatic state 1 (solid lines and solid squares) and state 
2 (dashed lines and open squares) as functions of the initial momentum Pq of the incident wave 
packet. 
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FIG. 3. Population difference as a function of fit with the following parameters: 6 = = 0.4, 
^ = 0.13, /3 = 12.5 and cjc = 1-0. (a) Comparison of results obtained with the FBTS and its 
variants with the exact quantum values, (b) Comparison of results obtained with the FBTS with 
proper and modified focused initial conditions. See text for definitions of the four cases. 
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FIG. 4. Population distribution of bacteriochlorophyll (Bchl.) 1-3 as function of t at temperature 
of 77 K. The solid lines represent the traceless-form results, and the corresponding color dots 
represent the trace- form results. The red data points are extracted from Ref. |34j 
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FIG. 5. Asymptotic adiabatic ground state population at 50 fs versus 7. Due to the strong 
system-bath couphng, the orthogonahty approximation compromises the FBTS results. 
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FIG. 6. The probability distribution Pg^s' for importance sampling of a pair of forward and back- 
ward quantum states versus t for a typical trajectory in the original basis (a) and rotated basis (b). 
In the later case, a more dispersed probability distribution allows for a more balanced samplings 
of all pairs of state combinations. One should replace 1/2 with +/- when using the legend in panel 
a to interpret the curves in panel b . 
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